/*
Project: SLMP impact evaluation 
Code to do the three-part balance tables (SLMPI, SLMPII and exterior control group).

*/
************************************************************************ 

clear all
		mata: mata clear
		cap log close
		set more off
		set trace off
		set linesize 140
		set memory 200m
		set matsize 1000
		set maxvar 120000
		set seed 123456789
		timer clear
		label drop _all 
		estimates clear
********************************************************************************
* DEFINITION OF THE WORKING DIRECTORY. INPUT AND OUTPUT FILES*

		if "`c(os)'"=="Windows" {								
		local drive: env HOMEDRIVE
		local path : env HOMEPATH
			global 	BASE_DIR "`drive'`path'\PATH_TO_MAIN_FOLDER" // change for path to main working directory
			}
		if "`c(os)'"=="MacOSX" {
		local path : env HOMEPATH
			global BASE_DIR "~//`path'\PATH_TO_MAIN_FOLDER" // change for path to main working directory
		}

*----------------------------
cd "$BASE_DIR"
 
 *IN*
		global INPUT	"$BASE_DIR\INPUT_DATA_PATH" //change for input data path

 *OUT*	
        global GRAPHS   "$BASE_DIR\OUTPUT_GRAPH_PATH"  //change for output graph path
	    global TABLES   "$BASE_DIR\OUTPUT_TABLES_PATH" //change for output tables path
		
****************************************************************************************************************

	use "$INPUT/SLMP_SIF_pixel_level", replace 

		drop if mws_pixel_id_5km==.

foreach var of varlist objectid mws_pixel_id_5km{
 tostring `var', replace 
 }
 
  
  	append using "$INPUT/Potential_controls_EAs_5km_pixel_level"

*------------------------------------------------------- 

 tostring year, gen(year_s) 
 tostring month, gen(month_s) 
 
	  gen season=.
	replace season=1 if (month>=10 | month==1)
	replace season=2 if  month>=2 & month<6
	replace season=3 if  month>=6 & month<10
 lab define season 1 "Dry Season" 2 "Pre-rainy season" 3 "Rainy season"

 la val season season 
  tostring season, gen(season_s)
  gen year_season_s=year_s+season_s
  destring year_season_s, gen(year_season)
  
    replace sif=. if sif<0.10
  preserve 
  sum sif
   gen nm=cond(sif!=.,0,1)
   bys mws_pixel_id_5km year_season: egen missing=total(nm)
   drop if missing!=0
   
 foreach var of varlist sif {
  bys mws_pixel_id_5km year_season: egen sif_s=mean(sif) // averaging at the year_season level
  }

  keep mws_pixel_id_5km year_season sif_s   
  bys mws_pixel_id_5km year_season: keep if _n==1
  tempfile average
  save `average', replace 
  restore 
   
  merge m:1 mws_pixel_id_5km year_season using `average', nogen
  

   keep if season==3 & sif_s!=. // just the pixels that effectively entered in the estimation 
   bys  mws_pixel_id_5km: keep if _n==1

  * replace  crop_mask_v2= crop_mask_02 if donor=="Control" & crop_mask_v2==.
   replace  rangeland_mask_v2 = range_mask_v2 if donor=="Control" & rangeland_mask_v2 ==.  
   replace  rangeland_mask_v2=0 if rangeland_mask_v2==. & donor=="Control"
   gen crop_range_mask = (crop_mask_v2+rangeland_mask_v2)
   
gen cat=.
replace cat=1 if donor=="WBI"
replace cat=2 if donor=="WBII"
replace cat=3 if donor=="Control"

global ylist "af250m_nut af250m_n_1 af250m_n_2  af_pH515cm af_ORCDRC1 af_EC515cm  af_Depthto af_CEC515c waterholdi elevation crop_mask_v2 rangeland_mask_v2 crop_range_mask crop_frac"

sum $ylist

g cat12=1 if cat==1
replace cat12=0 if cat==2

g cat13=1 if cat==1
replace cat13=0 if cat==3

g cat23=1 if cat==2
replace cat23=0 if cat==3

tokenize $ylist
while "`1'"!="" {
	sum `1' 
	local total=r(mean)
	
	reg `1' cat12 if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), robust cluster(mws_pixel_id_5km)
	local gr1=_b[cat12]+_b[_cons]
		local gr1=round(`gr1', 0.001)
	local gr2=_b[_cons]
		local gr2=round(`gr2', 0.001)
	local d12=_b[cat12]
		local d12=round(`d12', 0.001)

	mat v=e(V)
	local se12=sqrt(v[1,1])
		local se12=round(`se12', 0.001)
	local p12=normal(-abs(`d12')/`se12')*2
		local p12=round(`p12', 0.001)

	
	reg `1' cat13 if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), robust cluster(mws_pixel_id_5km)
	local gr3=_b[_cons]
		local gr3=round(`gr3', 0.001)
	local d13=_b[cat13]
		local d13=round(`d13', 0.001)
	mat v=e(V)
	local se13=sqrt(v[1,1])
		local se13=round(`se13', 0.001)
	local p13=normal(-abs(`d13')/`se13')*2
		local p13=round(`p13', 0.001)

	reg `1' cat23 if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), robust cluster(mws_pixel_id_5km)
	local d23=_b[cat23]
		local d23=round(`d23', 0.001)
	mat v=e(V)
	local se23=sqrt(v[1,1])
		local se23=round(`se23', 0.001)
	local p23=normal(-abs(`d23')/`se23')*2
		local p23=round(`p23', 0.001)

	oneway `1' cat if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), bon
	local dfm=r(df_m)
	local dfr=r(df_r)
	local fbon=r(F)
	local pbon=Ftail(`dfm', `dfr', `fbon')
		local pbon=round(`pbon', 0.001)
	
	mat output = nullmat(output)\((`gr1'), (`gr2'), (`gr3') ,(`d12') ,(`d13'),(`p12'), (`p13'))
	
	mac shift 
	}

mat colnames output =  group1 group2 group3 diff12 diff13 
mat rownames output = $ylist 
mat list output

count if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
local n=r(N)
count if cat==1 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
local n1=r(N)
count if cat==2 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
local n2=r(N)
count if cat==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
local n3=r(N)
mat output = (nullmat(output))\((`n1') ,(`n2'), (`n3'),(.) ,(.),(.) ,(.))
mat rownames output = $ylist N
mat list output


 	outtable using "$TABLES\Balance_three_groups", mat(output) caption("Balance analysis, three groups") replace


